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Abstract 

The two-dimensional random-bond Ising model is numerically studied on 
long strips by transfer-matrix methods. It is shown that the rate of decay of 
correlations at criticality, as derived from averages of the two largest Lyapunov 
exponents plus conformal invariance arguments, differs from that obtained 
through direct evaluation of correlation functions. The latter is found to 
be, within error bars, the same as in pure systems. Our results confirm 
field-theoretical predictions. The conformal anomaly c is calculated from the 
leading finite- width correction to the averaged free energy on strips. Estimates 
thus obtained are consistent with c = 1/2, the same as for the pure Ising 
model. 

PACS numbers: 05.50.-hq, 05.70.Jk, 64.60.Fr, 75.10.Nr 
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Conformal invariance methods have been important in the study of critical properties of 
two-dimensional spin models ^ , providing a number of exact results. Most applications have 
been to pure (translationally invariant) systems. When disorder is introduced, invariance 
under conformal transformations at criticality can be said to be preserved only in an average 
sense. This raises the question of whether averaged properties still connect to quantities in 
conformal theory, in the same way as their counterparts for pure systems. Further, it is well 
known that performing averages in inhomogeneous systems is a subtle matter 0-0, owing 
to the fact that the moments of local quantities scale with independent exponents. This is 
in contrast with the non-random case, where one has constant gap indices. 

In the present work, we study properties of the two-dimensional random-bond Ising 
model with a binary distribution of ferromagnetic interaction strengths, each occurring with 
equal probability. As the transition temperature is exactly known from duality one 
can be sure that numerical errors due to imprecise knowledge of the critical point are absent. 
In strip calculations, the only sources of such errors will then be the finite strip width and, 
for a random system, those arising from the averaging process. The former can be dealt 
with by finite-size scaling theory whereas one expects to reduce the effects of the latter 
by studying large enough samples. Averages of the following quantities are calculated at 
criticality: (i) the two largest Lyapunov exponents; [ii) spin-spin correlation functions; and 
(iii) free energy per spin. From (i) and conformal invariance arguments, we obtain what is 
supposed to represent the typical (as opposed to average) decay of the correlation 

functions. As shown below, results from (i) and (ii) differ in a way that corroborates 
the field-theoretic arguments of Ref. 0. Numerical data from {Hi) are used to estimate 



the leading finite-width correction to the bulk free energy, which is known pO| , pJ]] to be 
proportional to the conformal anomaly. 

We have used long strips of a square lattice, of width 2 < L < 13 sites with periodic 
boundary conditions. In order to provide samples that are sufficiently representative of 
disorder, we iterated the transfer matrix typically along 10^ lattice spacings. At each 
step, the respective vertical and horizontal bonds between first-neighbour spins i and j were 
drawn from a probability distribution 

1 

2' 

which ensures that the critical temperature j3c = l/ksTc of the corresponding two- 
dimensional system is given by 



P{Jij) = -Wr,-Jn)+5{J,,-rJ^)) , 0<r<l , (1) 



sinh(2/?eJo) sinh(2/5cr Jq) = 1 . (2) 

Two values of r were taken in calculations: r = 0.5 and 0.01. These roughly represent 
'weak' and 'strong' disorder, respectively. The critical temperatures, from Eq. 0, are: 
Tc (0.5)/ Jo = 1-641 . . .; (0.01)/ Jq = 0.5089 ... (to be compared with (l)/Jo = 2.269 . . .). 
The choice of r = 0.5 allowed us to compare our results with those of Glaus [|1^, who obtained 
averages of the two largest Lyapunov exponents and of the free energy for strips of L < 8. 
The reason for fixing r = 0.01 was that, in site-diluted Ising systems ||13[, a transfer-matrix 



ansatz shows that for concentrations such that Tc ~ 0.2Tc(l) one already is crossing over 
towards percolative behaviour. Analogy thus suggested that, for a critical temperature of 
the same order in our case, one might be able to detect similar effects (if they are present 
at all). 
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The procedure for evaluation of the largest Lyapunov exponent A^, for a strip of width 
L and length ^ 1 is well known |^,12]. Starting from an arbitrary initial vector Vq, one 
generates the transfer matrices % that connect columns i and i + 1, drawing bonds from the 
distribution Eq. [|, and applies them successively, to obtain: 

(3) 

The average free energy per site is, in units of /csT: 

fr\T) = -^Kl . (4) 

The second, third etc. exponents can generally be obtained by iteration of a set of initial 
vectors, Vj, orthogonal to each other and to vq; care must be taken to orthogonalise the 
iterated vectors every few steps 0], to prevent contamination with dominant-exponent com- 
ponents. Here, the task is easier since we shall be concerned only with the two largest 
Lyapunov exponents, and the Hamiltonian is invariant under global spin inversion. In order 
to calculate A° (A]^), it is then sufficient to iterate vq (vi) which is even (odd) under that 
same symmetry Hll2|, with no need for intermediate decontamination of the iterates of Vi. 



From A° and K\ one can form estimates of 

^-^ = ^l-^\ , (5) 

which, at the critical point, is related to the typical decay of correlations on a fixed sample 
in the plane ||^. The correlation functions in this case are predicted to have the form 

<aQaR>r~u R-^'^{A\yiR)-^I^ (fixed sample), (6) 

for ln(A In R) large, where A is proportional to the intensity of disorder. If now one tries to 
fit the standard relation between correlation length at criticality and the correlation-decay 
exponent ry: 

ri = L/niL , (7) 

which holds for pure systems, Eq. ^ implies that the result will in fact be an effective 
exponent, larger than 1/4 by an amount which increases with increasing disorder. As shown 
below, this is what arises from our numerical study. 

The direct calculation of correlation functions follows the lines of Section 1.4 of Ref. 
0, with standard adaptations for an inhomogeneous system. For two spins in, say, row 1, 
separated by a distance R, and for a given configuration C of bonds, one has: 

< >c= E ^(^o) ( n ^) / E ^(^o) ( n ^) ' (8) 

where (Tq = {ctq • • • } correspondingly for ct/j; the bonds that make the transfer matri- 
ces % belong to C . For pure systems the 2^— component vectors ip, ip are determined by the 
boundary conditions; for example, the choice of dominant left and right eigenvectors gives 
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the correlation function in an infinite system Here, one needs only be concerned with 
avoiding start-up effects, since there is no convergence of iterated vectors. This is done by 
discarding the first few hundred iterates of the initial vector vq. From then on, one can shift 
the dummy origin of Eq. ^ along the strip, taking -0 to be the left- (right-) iterate of 
Vq (vq ) up to the shifted origin. 

One performs averages of the correlation function (to be denoted by < ctqCt}^ > ) by 
shifting the origin along the strip and accumulating the corresponding results for Eq. |^. 
Such averages behave as 



< aoH, > ~ exp(-i?/e n . (9) 

For pure systems, ^ coincides with the definition oi in Eq. |^. Here we shall show that, 
at Tc one gets (within reasonable error bars) 

r''' = L/nv , (10) 

where = 1/4 is the same as for the pure system, with no logarithmic corrections such as 
those for typical correlations on a fixed sample (Eq. |^). This is again in agreement with the 
predictions of Ref . . 

The total strip length used was N = 10^ for L = 2 — 11 and 5 x 10^ for L = 12 
and 13 (except for free energy calculations, where N = 10^ was used also for the largest 
widths). Evaluation, both of Lyapunov exponents and of averaged correlation functions, 
involves iterating initial arbitrary vectors. In order to get rid of start-up effects, the first 
A^o = 2, 000 iterations were discarded. For Lyapunov exponents, the accumulated averages 
were recorded every 200 subsequent steps |12[. From this set we extracted estimates of 
exponents and corresponding fluctuations. This procedure was repeated three times with 
different initialisations. Final values of exponents and error bars were taken respectively as 
arithmetic averages of estimates at the end of each run (that is, each taking into account 
N — Nq steps), and root-mean-square averages of deviations (estimated from the sampling 
described above). We calculated correlation functions for distances R in Eq. || ranging from 
10 to 200 lattice spacings. The shifted origins were taken 200 columns apart, thus for each 
of the three independent runs, and each distance, one would have close to 500 samples of 
correlation functions (250 for L = 12 and 13). As all correlations are ferromagnetic, their 
fluctuations are well-behaved (compared e.g. to random-field problems, where frustration 
effects induce large-amplitude deviations from averages 



Our results are depicted in Figs. |l] (a) and (6), respectively for r = 0.5 and 0.01. We 
have plotted values of L/n^, both with ^ = given by Eq. ^ and ^ = ^ '^^^ of Eq. ^, against 
1/L^. The choice of horizontal axis is inspired by the pure case [15|, where it can be shown 
analytically that the leading finite-size corrections to rj are proportional to 1/L^. 

For weak disorder r = 0.5, Fig. ^ (a), both sequences are linear in 1/L^ to a good extent. 
Data from averaged correlation functions point inequivocally towards the pure-system value 
?7 = 1/4, in accordance with the predictions of Ref. [^], and also with very recent large-scale 
Monte Carlo simulations [0. Owing to the non-negligible fluctuations for large L, it is not 
feasible to try extrapolations of data for that region alone, as is done when no uncertainties 
are present |jl^. Instead, we turn to global fitting procedures. The analysis of least-squares 



fits of data (taking error bars into account) from L = Lq to L = 13, with Lq = 2, ... 10 



4 



indicates that stabilises at ~ 0.19 for Lq = 4 — 6; estimates of the extrapolated exponent 
and of errors are also very similar in this range. For these values of Lq a rough balance 
is attained, between the trend towards reducing errors by including fewer data, and that 
towards increasing uncertainties by giving more relative weight to fluctuations in individual 
data points. The exponent estimate from such fits is rj = 0.250 ± 0.002. 

Turning to the sequence of Lyapunov exponents for r = 0.5, an analysis of least-squares 
fits similar to that just described gives r] = 0.257 ± 0.002, with an even smaller ~ 0.05. 
If, following Ref. [0, one tries to plot data against 1/L, the amount of curvature is found 
to increase substantially. Taken on their own, the present results would seem to characterise 
a non- universal exponent, with leading finite-size corrections depending on When 
data for strong disorder are considered, however, this picture is not confirmed. As we shall 
see below, an overall scenario emerges in which the above estimate is to be regarded as an 
effective exponent, signalling the presence of logarithmic corrections of Eq. 

For strong disorder r = 0.01, Fig. |I| (6), the sequence of data from averaged correlation 
functions points again close to r/ = 1/4, though with a larger amount of scatter than in the 
previous case. Consideration of least-squares fits, analogously to that done above, leads to 
rj = 0.235 ± 0.015. The corresponding is ~ 1-8. 

Estimates of t] from Lyapunov exponents for r = 0.01 obviously do not scale linearly 
with 1/L^; one finds that the best linear fit is against L^*^ with (p ^ 0.2, in which case 
the extrapolate r]{L oo) is ~ 0. We refrain from attaching much significance to the 
latter result, as extrapolations with such low powers of 1/L are rather unreliable. Instead, 
we recall that, for a given finite L, estimates of rj are always larger than those for weak 
disorder. This is consistent with the logarithmic corrections for fixed-sample correlations 
displayed e.g. in Eq. Analytical expressions for the correlation length of energy-energy 
correlation functions p!B| , pi| , though not directly comparable with the present results, also 
show corrections in inverse powers of InL to the pure-system behaviour, whose absolute 
value increases with disorder. 

We now present results for the conformal anomaly. Table 1 shows the negative free energy 
per site, in units of ksT, as obtained from Eq. H, at the respective critical temperatures 
given by Eq. ^ for r = 1 (pure system), 0.5 and 0.01. Data for r = 1 are from Blote and 
Nightingale ||2^ , and shall be useful for comparison among different fitting procedures. The 
conformal anomaly c is related to the finite- width free energy at criticality by |]TD|,|TT|: 

TTC 

/l(TJ = /oo(T,) - — + ... . (11) 

We begin by recalling that, for pure Ising systems, finite-size estimates of c may approach 
the exact value c = 1/2 either from above or below, depending on what sort of corrections 
to the free energy are incorporated. If one stops at order 1/L^, the sequence of two-strip 
approximations with L and L + 1 reaches 1/2 from above; plugging in the next term (shown 
numerically by Blote and Nightingale pOf to be oc 1/L^), convergence of three-strip (L, L + 
1, L + 2) estimates is from below and much faster. This is illustrated in columns [i) and (n) 
of Table 2. The same happens when one performs straight-line or parabolic least-squares 
fits of data against l/L'^, using data from L = Lmm = 2, ... 11 up to L^ax = 13, as shown 
in columns {Hi) and [iv) of that table. 
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When disorder is considered, Table 1 shows that the error bars are a serious obstacle to 
evaluation of estimates from pairs or triplets of strip widths, especially for large L. In order 
to produce reliable results, one would need free-energy errors typically smaller than one part 
in 10^ for L > 10. The strips used would have to be respectively ~ 10^ (for r = 0.5) and 
10^ (for r = 0.01) times longer than the present value = 10^, assuming fluctuations to 
be normally distributed. Though not unfeasible, such task would demand a considerable 
amount of computer time. We shall then try to extract information from comparison between 
least-squares fits of disordered-system data and their pure-system counterparts. In Table 3 
are shown results of linear and parabolic fits against l/I/^, both for r = 0.5 and 0.01. In 
general, we found that fits starting at values of L larger than those displayed gave such large 
errors as to be rendered meaningless. For weak disorder, columns {i) and (ii), the overall 
trends closely resemble those found in the pure case of Table 2, including the approach 
of central estimates to c = 1/2 (from above or below, depending only on the type of fit), 
though of course with much larger errors. Results in column (i) are essentially the same as 



in Ref. |jT2|. As expected, error bars are even larger for strong disorder, columns {in) and 
(zf); though straight-line fits deteriorate very quickly, one still has reasonably well-behaved 
estimates from parabolas. However, if monotonic approach to c = 1/2 from below is present 
in the latter sequence, it is utterly smeared out by fiuctuations. 

Our results can be summarised as follows, 
(i) Averaged correlation functions at criticality of the disordered system decay with pure 
power-law behaviour, the same exponent = 1/4 of the pure Ising model. This is con- 
trary to early results, according to which disorder would lead to ?7 = PT[, but supports 



field-theoretical |^ and Monte-Carlo [|I^ evidence. Similarly, very recent results for site- 
diluted systems [^,^ indicate that, although susceptibility (7) and correlation length (z/) 
exponents seem to be concentration-dependent, the ratio / v = 2 — rj remains at 7/4 . 
{ii) Typical (as opposed to averaged) correlation decay, as obtained from the two largest Lya- 
punov exponents, behaves consistently with the presence of (non self-averaging) logarithmic 
corrections p. 

{Hi) Conformal anomaly estimates for the disordered Ising model behave in the same way 
(within error bars) as the corresponding ones for pure systems. We are thus led to state that 
no evidence has been found that 1/2. However, our results are not accurate enough to 
detect logarithmic corrections predicted by field theory |jl8|. 



One might view the estimate r] = 0.235 ± 0.015, from averaged correlation functions 
for r = 0.01, as signalling a crossover towards percolation-like behaviour (for which ||T5|j2^ 
rjp = 5/24). Though this is indeed expected |ll3| , |25| to occur sufficiently close to T = {i.e., 



r = 0), a systematic investigation of several points in that region is needed in order to sort 
out finite-size effects from r-dependent ones. We intend to do so, as part of a forthcoming, 
broader-ranging study |^ which will include generalisations of a calculational procedure 



initially devised for site-diluted problems [13| to random-bond and random-field systems 
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TABLES 



TABLE I. Critical free energies per site. Uncertainties in last quoted digits are shown in 
parentheses. 



T 

L 




r = 1 


2 


1 


886426125762 


3 


1 


842546256346 


4 


1 


828157728044 


5 


1 


821819028739 


6 


1 


818468940405 


7 


1 


816478784614 


8 


1 


815198160430 


9 


1 


814324882545 


10 


1 


813702482337 


11 


1 


813243149887 


12 


1 


812894445034 


13 


1 


812623456975 



U.o 






U.Ui 




1.0390 


(6) 


2 


121 


:4) 


0.9956 


(6) 


2 


081 


:4) 


0.9813 


(5) 


2 


068 


;4) 


0.9750 


(4) 


2 


062 


:3) 


0.9716 


(6) 


2 


059 


:4) 


0.9695 


(7) 


2 


056 


:5) 


0.9683 


(4) 


2 


055 


:3) 


0.9675 


(4) 


2 


055 


:3) 


0.9669 


(5) 


2 


054 


:4) 


0.9663 


(4) 


2 


053 


^3) 


0.9661 


(4) 


2 


053 


:3) 


0.9657 


(3) 


2 


052 


:3) 



TABLE II. Conformal anomaly estimates for pure Ising model, (i) : two-point fits (L, L + 1) 
with corrections, {ii) : three-point fits {L, L+l, L + 2) with and corrections. (Hi) : 
straight-line least-squares fits against with data from L to 13. (iv) : parabolic least-squares fits 
against with data from L to 13. Uncertainties in last quoted digits are shown in parentheses. 



L 


ii) 


(ii) 


(Hi) 




(iv) 




2 


0.60339151 


0.53003814 





5750 


(55) 


0.51113 


(211) 


3 


0.56530418 


0.49875466 





5408 


(39) 


0.49741 


(13) 


4 


0.53804550 


0.49507433 





5230 


(23) 


0.49722 


(27) 


5 


0.52348887 


0.49673455 





5146 


(14) 


0.49836 


(18) 


6 


0.51575527 


0.49823978 





5104 


(9) 


0.49909 


(10) 


7 


0.51133776 


0.49906829 





5080 


(6) 


0.49947 


(5) 


8 


0.50859245 


0.49948159 





5065 


(4) 


0.49966 


(3) 


9 


0.50676022 


0.49969077 





5055 


(3) 


0.49977 


(2) 


10 


0.50546906 


0.49980307 





5048 


(2) 


0.49983 


(1) 


11 


0.50452116 


0.49986769 





5042 


(2) 


0.49986769 


12 


0.50380296 







50380296 
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TABLE III. Conformal anomaly estimates for bond-disordered Ising model, (i) : straight-line 
least-squares fits against with data from L to 13; r = 0.5. {ii) : parabolic least-squares fits 

against 1/L^ with data from L to 13; r = 0.5. {Hi) : straight-line least-squares fits against 1/L^ 
with data from L to 13; r = 0.01. {iv) : parabolic least-squares fits against 1/L^ with data from 
L to 13; r = 0.01. Uncertainties in last quoted digits are shown in parentheses. 



L 


ii) 


{ii) 


{Hi) 


{iv) 


2 


0.569 (4) 


0.511 (3) 


0.53 (3) 


0.506 (10) 


3 


0.536 (10) 


0.498 (3) 


0.52 (7) 


0.508 (20) 


4 


0.521 (16) 


0.499 (6) 


0.51 (12) 


0.533 (45) 


5 


0.516 (23) 


0.495 (13) 


0.52 (17) 


0.503 (90) 


6 


0.509 (46) 


0.492 (26) 


0.51 (34) 




7 


0.502 (68) 









10 



FIGURES 

FIG. 1. rj = L/tt^, with ^ = of Eq. 5 (triangles) and ^ = ^ ^'"^ of Eq. 9 (circles) against 
The square on the vertical axis is at the pure system value 77 = 1/4 . Straight lines are 
least-square fits for L = 4 — 13. (a) : r = 0.5; {b) : r = 0.01 . 
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